Mo’orea Corallivores

Juvenile Chaetodon reticulatus. Image courtesy of Keoki Stender

Is there a relationship between decreasing coral cover and corallivore abundance in Backreef sites around Mo’orea, French Polynesia?

Coral reefs around the islands of Mo’orea are threatened by a variety of stressors such as cyclone activity and the outbreak of predatory sea stars (Acanthaster planci) that feed on coral tissue. As a response to these stressors and anthropogenic stressors, many locations around the island have undergone phase shifts in their benthic composition that involves the transformation of coral dominated communities to macroalgae dominated communities. In particular, the backreefs have been heavily transformed by the presence of Turbinaria, a well defended and persistent algae. I am interested if the resulting fish communities have changed as a result of this, particularly the abundance of Corallivores which feed primarily on scleractinian corals. In the following analysis, I will assess the impact of changing coral cover on corallivore populations.

The following analysis involves two datasets, both from the Mo’orea Coral Reef Long Term Ecological Research site.

fish : describe the species abundance and estimated size distributions (total body length to the greatest precision possible) of fishes surveyed as part of MCR LTER’s annual reef fish monitoring program. This dataset will be filtered to Backreef sites and fine trophic level of Corallivore.

coral_cover : describe the percentage cover of all stony corals (Scleractinia, pooled among genera) and other major groups analyzed from 0.5 x 0.5 m photographic quadrats at the Backreef habitat at the Moorea Coral Reef LTER, French Polynesia.

# Read in fish and coral data 
fish <- read_csv(here("data", "MCR_LTER_Annual_Fish_Survey_20250324.csv")) |> 
  janitor::clean_names()

coral_cover <- read_csv(here("data", "MCR_LTER_Coral_Cover_Backreef_Long_20250429.csv")) |> 
  janitor::clean_names()

Change in Corallivore Abundance through time

fish_2006 <- fish |> 
  group_by(fine_trophic, habitat) |> 
  filter(fine_trophic != "na")|> 
  filter(year == "2006") |> 
  filter(habitat == "Backreef") |> 
  summarise(n = n()) |> 
  mutate(total_number = n) |> 
  arrange(desc(total_number))
          
          
fish_2006_plot <- ggplot(fish_2006, aes(x = total_number, 
                      y = reorder(fine_trophic, +total_number), 
                      fill = fine_trophic)) +
  geom_col() + 
  ylab(" ") + 
  xlab("Count") + 
  labs(title = " 2006 Backreef Fish Abundance by Group", 
       x = "Count") +
  theme(legend.position = "none")

fish_2024 <- fish |> 
  group_by(fine_trophic, habitat) |> 
  filter(fine_trophic != "na")|> 
  filter(year == "2024") |> 
  filter(habitat == "Backreef") |> 
  summarise(n = n()) |> 
  mutate(total_number = n) |> 
  arrange(desc(total_number))

fish_2024_plot <- ggplot(fish_2024, aes(x = total_number, 
                      y = reorder(fine_trophic, +total_number), 
                      fill = fine_trophic)) +
  geom_col() + 
  ylab(" ") + 
  xlab("Count") + 
  labs(title = " 2024 Backreef Fish Abundance by Group", 
       x = "Count") +
  theme(legend.position = "none")

Corallivores fell from the third‑most abundant group in the first survey year (202 individuals) to the sixth‑most abundant group in the most recent sampling year (79 individuals) at the LTER Backreef sites. The following figure shows the drop from 2006 to 2024.

corallivores_through_years <- fish |> 
  group_by(year) |> 
  filter(habitat == "Backreef") |> 
  filter(fine_trophic == "Corallivore") |> 
  summarise(n = n()) |> 
  mutate(total_number = n)
  
corallivore_plot <- ggplot(corallivores_through_years, 
       aes(x = year, 
           y = total_number)) + 
  geom_line(color = "lightblue") + 
  labs(x = "Year", 
    y = "Total Number of Fish", 
    title = "Change in Backreef Corallivore Abundance (2006-2024)")

corallivore_plot

Relationships and causal relationships described with a DAG

The three predictor variables that will be included in this model are LTER site, year, and mean coral cover of transects within a site. The site and year should determine the coral cover which in turn should determine the response; proportion corallivore in the fish community.

Statistical model

\[ \begin{align} \text{Proportion Corallivores} &\sim Beta(\mu, \phi) \\ logit(\mu) &= \beta_0 + \beta_1 \text{Mean Site Coral Cover} + \beta_2 \text{Site} + \beta_3 \text{Year} \end{align} \] The standard parameters of beta distribution are α and β, but for beta regression, it is more useful to reparameterize them to μ and ϕ.

For future reference:

α = μ*ϕ

β = (1- μ)*ϕ

The response variable is a proportion between 0 and 1. We want to investigate how that changes with the changes in the predictor variables.

μ: the mean of beta distribution.

ϕ: precision parameter. The larger the ϕ, less variance in distribution. The inverse is true: a low ϕ, yields more variance in distribution.

Because μ needs to be between 0 and 1, a transformation using the logit link is necessary to convert values to that range.

Simulate data according to model assumptions / Model fit to simulated data recovers the parameters

The following parameters will be assigned these values:

β0 = .5

β1 = 1

ϕ = 5

We should recover them after the following data simulation and model fit.

# Sample size for simulation  
n <- 5000
# Precision parameter
phi <- 5
# Intercept 
beta0 <- .5 
# Predictor coefficient 
beta1 <- 1 


# Generate random numbers from a normal distribution of sample size 
predictor <- rnorm(n)

# Logit mu regression line 
logit_mu <- beta0 + beta1 * predictor 

# plogis() converts logit_mu to mu which will be bound between 0 and 1
mu <- plogis(logit_mu) 

# Alpha shape of beta distribution or shape 1 
alpha <- mu * phi
# Beta shape of beta distribution or shape 2 
beta  <- (1 - mu) * phi


# Simulate response variable with sample size n and two shape parameters using rbeta()
response <- rbeta(n, alpha, beta)

# Fit beta regression model to simulated data 
model <- betareg(response ~ predictor)

# Obatain parameters: b0, b1, and phi 
summary(model)

Call:
betareg(formula = response ~ predictor)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-3.4755 -0.6639  0.0159  0.6858  4.0502 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.49657    0.01237   40.15   <2e-16 ***
predictor    1.00919    0.01370   73.64   <2e-16 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  4.97625    0.09528   52.23   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood:  2611 on 3 Df
Pseudo R-squared: 0.5095
Number of iterations: 12 (BFGS) + 2 (Fisher scoring) 

The summary returned the following:

Intercept (β0): 0.49657

predictor (β1): 1.00919

phi(ϕ): 4.97625

All values were significant and extremely close to the originally assigned numbers above.

Inference

Alternative Hypothesis: There is a positive relationship between mean site coral cover and proportion corallivore in the backreef fish community.

Null Hypothesis : There is no relationship between mean site coral cover and proportion corallivore in the backreef fish community.

First, processing data to create a dataframe which contains the corallivore proportion at each LTER backreef site each year, and the mean coral cover of monitored transects at each LTER backreef site, each year.

corallivore <- fish |> 
  filter(habitat == "Backreef") |> 
  group_by(year, site ) |> 
  summarise(proportion_corallivore = mean(fine_trophic == "Corallivore", na.rm = TRUE)) 

coral_cover <- coral_cover |> 
  filter(benthic_category == "coral") |> 
  group_by(year, site) |> 
  summarise(mean_site_coral_cover = mean(percent_cover)) 

full_data_frame <- left_join(coral_cover, corallivore, join_by(year, site)) 

full_data_frame <- full_data_frame |> 
  filter(year != 2005)

Plot relationship

site_labels <- c(
  "LTER_1" = "LTER 1",
  "LTER_2" = "LTER 2",
  "LTER_3" = "LTER 3",
  "LTER_4" = "LTER 4",
  "LTER_5" = "LTER 5",
  "LTER_6" = "LTER 6"
)

ggplot(full_data_frame, 
       aes(x = mean_site_coral_cover, 
           y = proportion_corallivore, 
           color = year)) + 
  scale_color_gradient(low = "lightblue", high = "darkblue") + 
  geom_point() + 
  facet_wrap("site",labeller = labeller(site = site_labels)) + 
  labs(x = "Mean Site Coral Cover (%)", 
       y = "Proportion Corallivore", 
       color = "Year")

# Fit beta regression model to data using glmmTMB()

corallivore_model <- glmmTMB(proportion_corallivore ~ mean_site_coral_cover + site + year,
  data = full_data_frame,
  family = beta_family(link = "logit"))

summary(corallivore_model)

 Family: beta  ( logit )
Formula:          proportion_corallivore ~ mean_site_coral_cover + site + year
Data: full_data_frame

      AIC       BIC    logLik -2*log(L)  df.resid 
   -560.1    -535.5     289.1    -578.1       105 


Dispersion parameter for beta family ():  175 

Conditional model:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)           75.265502  11.098162   6.782 1.19e-11 ***
mean_site_coral_cover  0.008127   0.003320   2.448  0.01435 *  
siteLTER_2             0.344175   0.109374   3.147  0.00165 ** 
siteLTER_3             0.865243   0.102320   8.456  < 2e-16 ***
siteLTER_4             0.593038   0.119589   4.959 7.09e-07 ***
siteLTER_5             0.595067   0.104723   5.682 1.33e-08 ***
siteLTER_6             0.653496   0.105820   6.176 6.59e-10 ***
year                  -0.038947   0.005493  -7.090 1.34e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Create grid of with percentage values of coral cover from 0 to 1. Reference year is 2006, reference site is LTER 1

pred_grid <- expand_grid(mean_site_coral_cover = seq(0,1, by = 0.01)) |> 
  mutate(site = "LTER_1", 
    year  = 2006)

# Make predictions with corallivore model along with standard error of each response value

pred_se <- predict(
  corallivore_model,
  newdata = pred_grid,
  type = "response",
  se.fit = TRUE
)

# Combine predictions with upper/lower bounds of confidence intervals

prediction_final <- pred_grid |> 
  mutate(
    fit = pred_se$fit,
    se = pred_se$se.fit,
    lower = fit - 1.96 * se,
    upper = fit + 1.96 * se
  )

# Plot expected response of model and CIs

ggplot(prediction_final, aes(x = mean_site_coral_cover, y = fit)) +
  geom_line(colour = "lightblue") +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              fill = "lightblue", alpha = 0.2) +
  scale_y_continuous(limits = c(0, 0.10),   # y‑axis from 0 to 0.10
                     expand = expansion(mult = c(0, 0.05))) +
  labs(
    x = "Mean site coral cover",
    y = "Proportion corallivore"
  ) +
  theme_minimal()

A hypothesis is tested and the evidence is interpreted

The coefficients obtained from the beta regression model represent the change in the log-proportion of the response variable for a one-unit change in the predictor variable. In other words, for every one percent increase in mean site coral cover, we would expect the proportion of corallivore in the fish community to increase by 0.008127. The coefficients for mean site coral cover, year , and all sites were significant. We can successfully reject the null hypothesis. Our data is consistent with the alternative hypothesis that there is a positive relationship between mean site coral cover and proportion corallivore in the backreef fish community.